  # This program is AIM/Enduse and CGE Japan coupling analysis
  
  library(gdxrrw)
  library(ggplot2)
  library(dplyr)
  library(reshape2)
  library(tidyr)
  library(maps)
  library(grid)
  library(RColorBrewer)
  library(gridExtra)
  library(tableone)
  OrRdPal <- brewer.pal(9, "OrRd")
  YlGnBupal <- brewer.pal(9, "YlGnBu")
  Redspal <- brewer.pal(9, "Reds")
  pastelpal <- brewer.pal(9, "Pastel1")
  pastelpal <- brewer.pal(8, "Set1")
  dark2pal <- brewer.pal(8, "Dark2")
  
  MyThemeLine <- theme_bw() +
    theme(
      panel.border=element_rect(fill=NA),
      panel.grid.minor = element_line(color = NA), 
      axis.line=element_line(colour="black"),
      panel.background=element_rect(fill = "white"),
      panel.grid.major=element_blank(),
      strip.background=element_rect(fill="white", colour="white"),
      strip.text.x = element_text(size=12, colour = "black", angle = 0,face="bold"),
      axis.text.x=element_text(angle=45, vjust=0.9, hjust=1, margin = unit(c(t = 0.3, r = 0, b = 0, l = 0), "cm")),
      axis.text.y=element_text(margin = unit(c(t = 0, r = 0.3, b = 0, l = 0), "cm")),
      legend.text = element_text(size = 7),
      legend.title = element_text(size = 7),
      axis.ticks.length=unit(-0.15,"cm"),
      strip.text.y = element_text(size = 12)
    )
  
  
  #-- data load
  dir.create("../output/")
  outputdir <- c("../output/")
  linepalette <- c("#4DAF4A","#FF7F00","#377EB8","#E41A1C","#984EA3","#FFFF33","#A65628","#F781BF","#8DD3C7","#FB8072","#80B1D3","#FDB462","#B3DE69","#FCCDE5","#D9D9D9","#BC80BD","#CCEBC5","#FFED6F","#7f878f")
  landusepalette <- c("#8DD3C7","#FF7F00","#377EB8","#4DAF4A","#A65628")
  scenariomap <- read.table("../data/scenariomap.map", sep="\t",header=T, stringsAsFactors=F)
  Endusescenarios <- c("Enduse1","Enduse2","Enduse3","Enduse4","Enduse5","Enduse2E0p5","Enduse2E2p0")
  CGEload0 <- rgdx.param('../data/JPN_emf.gdx','EMFtemp1') %>% rename("Value"=EMFtemp1,"Variable"=VEMF) %>% mutate(MODEL="CGE")
  EUload0.0 <- rgdx.param('../data/AIMEnduse.gdx','EMFtemp1') %>% rename("SCENARIO"=i1,"REMF"=i2,"Variable"=i3,"Y"=i4,"Value"=value) %>% mutate(MODEL="Enduse")
  
  EUload0 <- spread(EUload0.0,key="SCENARIO",value=Value)
  EUload1 <- gather(EUload0,key="SCENARIO",value="Value",-REMF,-Variable,-Y,-MODEL)
  allmodel.1 <- inner_join(rbind(CGEload0,EUload1),scenariomap,by="SCENARIO")
  allmodel.1$Y <- as.numeric(levels(allmodel.1$Y))[allmodel.1$Y]
  
  #-- Policy cost
  PolCost0 <- allmodel.1 %>% filter(Model %in% Endusescenarios & Variable %in% c("Pol_Cos_Add_Tot_Ene_Sys_Cos","GDP_MER"))
  PolCost1 <- spread(PolCost0,key=Variable,value=Value)                                
  PolCost1$Pol_Cos_GDP_Los_rat <- PolCost1$Pol_Cos_Add_Tot_Ene_Sys_Cos/PolCost1$GDP_MER*100
  PolCost2 <- gather(PolCost1,key=Variable,value=Value,-REMF,-SCENARIO,-Y,-MODEL,-Model,-CLP) %>% filter(Variable=="Pol_Cos_GDP_Los_rat") 
  allmodel.2 <- rbind(allmodel.1,PolCost2)
  allmodel <- allmodel.2 %>% filter(Model %in% c("CGE1","CGE2","Enduse1"))
  #Prepareation of object for power point
  varlist <- read.table("../data/varlist.txt", sep="\t",header=F, stringsAsFactors=F)
  
  arealist <- data.frame(c("TPES","POWER","TFC"),c("EJ/year","EJ/year","EJ/year"),c("Primary energy","Power generation","Final energy"))
  names(arealist) <- names(varlist)
  nalist <- rbind(varlist,arealist)
  allplot <- as.list(nalist)
  plotflag <- as.list(nalist)
  data4plot <- as.list(nalist)

  #for (i in 1:nrow(varlist)){
  for (i in 1:nrow(varlist)){
    if(length(allmodel[allmodel$Variable==varlist[i,1],c(1)])>0){
      data4plot[[i]] <- allmodel %>% filter (Variable==varlist[i,1] & Model %in% c("CGE1","CGE2","Enduse1")) 
      plot.0 <- ggplot() + 
        geom_line(data=data4plot[[i]][data4plot[[i]]$Variable==varlist[i,1]  & data4plot[[i]]$CLP %in% c("Baseline","Mitigation"),c(-2)],aes(x=Y, y = Value , color=Model,group=interaction(Model)),stat="identity") +
        geom_point(data=data4plot[[i]][data4plot[[i]]$Variable==varlist[i,1] & data4plot[[i]]$CLP %in% c("Baseline","Mitigation"),c(-2)],aes(x=Y, y = Value , color=Model,group=interaction(Model),shape=Model),size=3.0,fill="white") +
        MyThemeLine+ scale_color_manual(values=linepalette,labels = c(CGE1="CGE Stand-alone",CGE2="Integrated model",Enduse1="Enduse Stand-alone")) + 
        scale_shape_manual(values=c(1,2,3,23),labels = c(CGE1="CGE Stand-alone",CGE2="Integrated model",Enduse1="Enduse Stand-alone")) +
        xlab("year") + ylab(varlist[i,3])  +  ggtitle(varlist[i,2]) +
      annotate("segment",x=2005,xend=2050,y=0,yend=0,linetype="dashed",color="grey")+ 
        theme(legend.title=element_blank()) 
      plot.1 <- plot.0 + facet_wrap(~ CLP,scales="free")
      outname <- paste0(outputdir,"",varlist[i,1],".png")
      ggsave(plot.1, file=outname, dpi = 150, width=7, height=4,limitsize=FALSE)
    allplot[[i]] <- plot.1
    }
    plotflag[[i]] <- nrow(allmodel[allmodel$Variable==varlist[i,1] ,c(1)])
  }
  
  tpespalette <- c("Coal|w/o CCS"="#000000","Coal|w/ CCS"="#7f878f","Oil|w/o CCS"="#ff2800","Oil|w/ CCS"="#ffd1d1","Gas|w/o CCS"="#9a0079","Gas|w/ CCS"="#c7b2de","Hydro"="#0041ff","Nuclear"="#663300","Solar"="#b4ebfa","Wind"="#ff9900","Biomass|w/o CCS"="#35a16b","Biomass|w/ CCS"="#cbf266","Geothermal"="#edc58f","Other"="#ffff99",
                   "Solid|Biomass"=dark2pal[1],"Solid|Coal"=dark2pal[2],"Liquid|Biomass"=dark2pal[3],"Liquid|Oil"=dark2pal[4],"Gas"=dark2pal[5],"Electricity"=dark2pal[6],"Heat"=dark2pal[7])
#  "Solid|Biomass"="#cbf266","Solid|Coal"=pastelpal[1],"Liquid|Biomass"="#35a16b","Liquid|Oil"=pastelpal[2],"Gas"=pastelpal[3],"Electricity"=pastelpal[4],"Heat"=pastelpal[5])
  
  areamap <- read.table("../data/Areafigureorder.txt", sep="\t",header=T, stringsAsFactors=F)
  areamappara <- read.table("../data/Area.map", sep="\t",header=T, stringsAsFactors=F)
  area.0 <- allmodel %>% filter(Variable %in% as.vector(areamap$Variable)) %>% left_join(areamap,by="Variable") %>% ungroup()
  labeli <- as_labeller(c(`CGE2` = "Integrated model",
                          `CGE1` = "CGE stand-alone",
                          `Enduse1` = "Enduse Stand-alone",
                          'Baseline'="Baseline",
                          'Mitigation'="Mitigation"))
  plot.1 <- function(){
    plot <- ggplot() + geom_area(data=XX[XX$CLP %in% c("Baseline","Mitigation"),c(-100)],aes(x=Y, y = Value , fill=reorder(Ind,-order)), stat="identity") + 
      ylab(ylab1) + xlab(xlab1) +labs(fill="")+ guides(fill=guide_legend(reverse=TRUE)) + MyThemeLine +
      theme(legend.position="bottom", text=element_text(size=12),  
            axis.text.x=element_text(angle=0, vjust=0.9, hjust=1, size = 12)) +
      guides(fill=guide_legend(ncol=5))
    plot2 <- plot + facet_grid(CLP  ~ Model, labeller=labeli,scales="free") + scale_fill_manual(values=colorpal)+
      annotate("segment",x=2005,xend=2050,y=0,yend=0,linetype="solid",color="grey") + theme(legend.position='bottom')
    return(plot2)
  }
  
  for(j in 1:length(areamappara)){
    XX <- area.0 %>% filter(Class==areamappara[j,1]) %>% select(Model,CLP,Ind,Y,Value,order)  %>% arrange(order)
    na.omit(XX$Value)
    data4plot[[i+j]] <- XX
    unit_name <-areamappara[j,3]
    ylab1 <- paste0(areamappara[j,2], " (", unit_name, ")")
    xlab1 <- areamappara[j,2]
    colorpal <- tpespalette
    plot_TPES.1 <- plot.1()
    allplot[[i+j]] <- plot_TPES.1
    outname <- paste0(outputdir,areamappara[j,1],".png")
    ggsave(plot_TPES.1, file=outname, dpi = 450, width=9, height=6,limitsize=FALSE)
    
  }
  
  #Decomposition analysis
  for(i in 1:2){
    if(i==1){
      Decom1 <- rgdx.param('../data/analysis.gdx','Loss_dcp_gdp' ) %>% rename("value"=Loss_dcp_gdp,"Sector"=SCO2_S,"Element"=decele) %>% select(-R)
    }else{
        Decom1 <- rgdx.param('../data/analysis.gdx','Loss_dcp_rate' ) %>% rename("value"=Loss_dcp_rate,"Sector"=SCO2_S,"Element"=decele) %>% select(-R)
      Decom1 <- Decom1 %>% filter(!Sector %in% c("BIO","CCS","OEN","FFE","PWR"))
    }
    Decom2 <- Decom1 %>% filter(SCENARIO %in% c("SSP2_AIMEInv2_80NPi_NoCC","SSP2_80NPi_NoCC") & Y %in% c(2030,2050) & Sector!="OTH" & Sector!="tot") 
  levels(Decom2$SCENARIO)[levels(Decom2$SCENARIO)=="SSP2_AIMEInv2_80NPi_NoCC"] <- "Integrated model"
  levels(Decom2$SCENARIO)[levels(Decom2$SCENARIO)=="SSP2_80NPi_NoCC"] <- "CGE stand-alone"
  levels(Decom2$Element)[levels(Decom2$Element)=="del_output"] <- "Output change"
  levels(Decom2$Element)[levels(Decom2$Element)=="outputratio"] <- "Value-added_Output ratio"
  levels(Decom2$Element)[levels(Decom2$Element)=="residual"] <- "Residual"
  flabel <- c("Change in % of GDP","sectors")
    plot <- ggplot() + geom_bar(data=Decom2[Decom2$Element!="va",c(-20)],aes(x=Sector, y = value*100 , fill=Element), stat="identity") +
      geom_point(data=Decom2[Decom2$Element=="va",c(-20)],aes(x=Sector, y = value*100 ),color="black", stat="identity") +
      ylab(flabel[1]) + xlab(flabel[2]) +labs(fill="") + scale_fill_manual(values=landusepalette) +
      guides(fill=guide_legend(reverse=TRUE)) + 
      MyThemeLine + 
      theme(legend.position="bottom", text=element_text(size=12),legend.text=element_text(size=12)) +
      guides(fill=guide_legend(ncol=5))
    plot2 <- plot + facet_grid(SCENARIO~Y,scales="free_x") + 
      annotate("segment",x=0,xend=10,y=0,yend=0,linetype="dashed",color="grey") 
      
    outname <- paste0(outputdir,"decomp",i,".png")
    ggsave(plot2, file=outname, dpi = 450, width=7, height=5,limitsize=FALSE)
  
  }
  
#--- Degree of coincidence
  system(paste("gams dataprocess.gms",sep=" "))

#---- For paper figure 
#
plot_fig1 <- list(1:10)
j <- 1
techist <- c("","_NoNuc","_NoCCS")
ilist <- c(74,77)
ColorVecSce <- c(linepalette[1],linepalette[2])
ShapeVecSce <- c(1,2)
layout2 <- rbind(c(1, 2, 3,4),
                 c(1, 2, 3,5))

for(l in 1:length(techist)){
  dumrow <- data.frame("SSP2_80NPi_NoCC","JPN","Prc_Car","2000","0","CGE","CGE2",paste0("Baseline",techist[l]))
  names(dumrow) <- names(data4plot[[77]])
  dumrow$Y <- as.numeric(levels(dumrow$Y))[dumrow$Y]
  dumrow$Value <- as.numeric(levels(dumrow$Value))[dumrow$Value]
  data4plot[[77]] <- rbind(data4plot[[77]],dumrow)
  for(j in 1:3){
    unit_name <-areamappara[j,3]
    ylab1 <- paste0(areamappara[j,2], " (", unit_name, ")")
    xlab1 <- areamappara[j,2]
    plotdata1 <- data4plot[[nrow(varlist)+j]] %>% filter(Model %in% c("CGE2") & CLP %in% c(paste0("Baseline",techist[l]),paste0("Mitigation",techist[l])))
    plot <- ggplot() + geom_area(data=plotdata1,aes(x=Y, y = Value , fill=reorder(Ind,-order)), stat="identity") + 
      ylab(ylab1) + xlab("") +labs(fill="")+ MyThemeLine +
      theme(legend.position="bottom", text=element_text(size=12),  
            axis.text.x=element_text(vjust=0.9, hjust=1, size = 12)) +
      guides(fill=FALSE)
      #    guides(fill=guide_legend(ncol=5)) + guides(fill=guide_legend(reverse=TRUE)) 
    plot_fig1[[j]] <- plot + facet_grid(CLP ~ ., scales="free") + scale_fill_manual(values=colorpal)+
      annotate("segment",x=2005,xend=2050,y=0,yend=0,linetype="solid",color="grey") + 
      theme(strip.text.x = element_blank(),strip.text.y = element_blank()) # Remove label
  }

  
  for(k in 1:length(ilist)){
    i <- ilist[k]
    dataplot1 <- data4plot[[i]][data4plot[[i]]$Variable==varlist[i,1] & data4plot[[i]]$Model=="CGE2",c(-2)] %>% filter(CLP %in% c(paste0("Baseline",techist[l]),paste0("Mitigation",techist[l])))
    plot_fig1[[j+k]] <- 
      ggplot() + 
    geom_line(data=dataplot1,aes(x=Y, y = Value , color=CLP,group=interaction(CLP)),stat="identity") +
    geom_point(data=dataplot1,aes(x=Y, y = Value , color=CLP,group=interaction(CLP),shape=CLP),size=3.0,fill="white") +
    MyThemeLine+ scale_color_manual(values=ColorVecSce) + scale_shape_manual(values=ShapeVecSce) +
    xlab("") + ylab(varlist[i,3])  + xlim(2005,2050) +
    annotate("segment",x=2005,xend=2050,y=0,yend=0,linetype="dashed",color="grey")+ 
    theme(legend.position="bottom", text=element_text(size=12),  
            axis.text.x=element_text(vjust=0.9, hjust=1, size = 12)) +
      guides(color=FALSE,shape=FALSE)
    #    theme(legend.title=element_blank()) 
    if (k ==1){
      plot_fig1[[j+k]] <- plot_fig1[[j+k]]+
        theme(axis.text.x=element_blank()) 
    }
  }
  outname <- paste0("../output/4paper/Fig1",techist[l],".png")
  png(filename=outname, width=4000, height=2000,res=500)
  grid.arrange(plot_fig1[[1]], plot_fig1[[2]], plot_fig1[[3]], plot_fig1[[4]],plot_fig1[[5]],
               layout_matrix = layout2)
  dev.off()  
}
  
#GDP loss rate  
  plot_fig4 <- as.list(5)

  gdplossname <- varlist %>% filter(V1=="Pol_Cos_GDP_Los_rat")
  plotdata1 <- allmodel %>% filter (Variable==gdplossname[,1])
  plotdata1 <- na.omit(plotdata1)
  
  plotdata0 <- plotdata1 %>% filter(Model %in% c("CGE1","CGE2","Enduse1") & CLP==c("Mitigation")) %>% select(-SCENARIO)
  plot_fig4[[1]] <- 
    ggplot() + 
    geom_line(data=plotdata0,aes(x=Y, y = Value , color=Model,group=interaction(Model)),size=1.0) +
    geom_point(data=plotdata0,aes(x=Y, y = Value , color=Model,group=interaction(Model),shape=Model),size=3.0,fill="white") +
    MyThemeLine+ scale_color_manual(values=linepalette,labels = c(CGE1="CGE Stand-alone",CGE2="Integrated model",Enduse1="Enduse Stand-alone")) +
    scale_shape_manual(values=c(1,2,3,23),labels = c(CGE1="CGE Stand-alone",CGE2="Integrated model",Enduse1="Enduse Stand-alone")) +
#    scale_shape_manual(values=ShapeVecSce) +
    xlab("Year") + ylab(paste0("GDP loss rates (",gdplossname[,3],")") ) + ylim(-0.1,3.5) +
    annotate("segment",x=2005,xend=2050,y=0,yend=0,linetype="dashed",color="grey")+ 
    theme(legend.position="bottom", text=element_text(size=12), legend.title=element_blank(),
          axis.text.x=element_text(vjust=0.9, hjust=1, size = 12)) +
    guides(shape=guide_legend(ncol=1),color=guide_legend(ncol=1)) +
    ggtitle("a        ")
  
  plotdata1_tech0 <- plotdata1 %>% filter(Model %in% c("CGE1","CGE2")) %>% select(-SCENARIO) 
  plotdata1_tech0 <- na.omit(plotdata1_tech0)
  plotdata1_tech <- spread(plotdata1_tech0,key=Model,value=Value)
  plot_fig4[[2]] <- 
    ggplot() + 
    geom_line(data=plotdata1_tech,aes(x = CGE1, y = CGE2 , color=CLP),size=1) +
    geom_point(data=plotdata1_tech,aes(x = CGE1, y = CGE2 , color=CLP,shape=CLP),size=3.0,fill="white") +
    MyThemeLine+ scale_color_manual(values=pastelpal,labels = c(Mitigation="Default",Mitigation_NoCCS="No CCS",Mitigation_NoNuc="No nuclear")) +
    scale_fill_manual(values=pastelpal,labels = c(Mitigation="Default",Mitigation_NoCCS="No CCS",Mitigation_NoNuc="No nuclear")) +
    scale_shape_manual(values=c(17,18,19),labels = c(Mitigation="Default",Mitigation_NoCCS="No CCS",Mitigation_NoNuc="No nuclear")) +
    ylab(paste0("Integrated model (",gdplossname[,3],")")) + xlab(paste0("CGE stand-alone (",gdplossname[,3],")") )  +
    annotate("segment",x=0,xend=4,y=0,yend=4,linetype="dashed",color="black")+ 
    theme(legend.position="bottom", text=element_text(size=12), legend.title=element_blank(),
          axis.text.x=element_text(vjust=0.9, hjust=1, size = 12)) +
    guides(shape=guide_legend(ncol=1),color=guide_legend(ncol=1)) +
    ggtitle("c     ")
  
  plotdata20 <- plotdata1 %>% filter(Model %in% c("Enduse1","CGE2")) %>% select(-SCENARIO,-MODEL) 
  plotdata20 <- na.omit(plotdata20)
  plotdata2 <- spread(plotdata20,key=Model,value=Value)
  plot_fig4[[3]] <- 
    ggplot() + 
    geom_line(data=plotdata2,aes(x = Enduse1, y = CGE2 , color=CLP),size=1.0) +
    geom_point(data=plotdata2,aes(x = Enduse1, y = CGE2 , color=CLP,shape=CLP,fill=CLP),size=3.0) +
    MyThemeLine+ 
    scale_color_manual(values=pastelpal,labels = c(Mitigation="Default",Mitigation_NoCCS="No CCS",Mitigation_NoNuc="No nuclear")) +
    scale_fill_manual(values=pastelpal,labels = c(Mitigation="Default",Mitigation_NoCCS="No CCS",Mitigation_NoNuc="No nuclear")) +
    scale_shape_manual(values=c(17,18,19),labels = c(Mitigation="Default",Mitigation_NoCCS="No CCS",Mitigation_NoNuc="No nuclear")) +
    ylab(paste0("Integrated model (",gdplossname[,3],")")) + xlab(paste0("Energy system model (",gdplossname[,3],")") ) +  
    annotate("segment",x=0,xend=1.5,y=0,yend=1.5,linetype="dashed",color="black")+ 
    theme(legend.position="bottom", text=element_text(size=12), legend.title=element_blank(),
          axis.text.x=element_text(vjust=0.9, hjust=1, size = 12)) +
    guides(shape=guide_legend(ncol=1),color=guide_legend(ncol=1)) +
    ggtitle("d      ")

# Equivalent variation
  gdplossname2 <- varlist %>% filter(V1=="Pol_Cos_Equ_Var_rat")
  plotdata1 <- allmodel %>% filter (Variable==gdplossname2[,1])
  plotdata1 <- na.omit(plotdata1)
  
  plotdata0 <- plotdata1 %>% filter(Model %in% c("CGE1","CGE2","Enduse1") & CLP==c("Mitigation")) %>% select(-SCENARIO)
  plot_fig4[[4]] <- 
    ggplot() + 
    geom_line(data=plotdata0,aes(x=Y, y = Value , color=Model,group=interaction(Model)),size=1.0) +
    geom_point(data=plotdata0,aes(x=Y, y = Value , color=Model,group=interaction(Model),shape=Model),size=3.0,fill="white") +
    MyThemeLine+ scale_color_manual(values=linepalette,labels = c(CGE1="CGE Stand-alone",CGE2="Integrated model",Enduse1="Enduse Stand-alone")) +
    scale_shape_manual(values=c(1,2,3,23),labels = c(CGE1="CGE Stand-alone",CGE2="Integrated model",Enduse1="Enduse Stand-alone")) +
    #    scale_shape_manual(values=ShapeVecSce) +
    xlab("Year") + ylab(paste0("Equivalent variation (",gdplossname2[,3],")") ) + ylim(-0.1,3.5) +
    annotate("segment",x=2005,xend=2050,y=0,yend=0,linetype="dashed",color="grey")+ 
    theme(legend.position="bottom", text=element_text(size=12), legend.title=element_blank(),
          axis.text.x=element_text(vjust=0.9, hjust=1, size = 12)) +
    guides(shape=guide_legend(ncol=1),color=guide_legend(ncol=1)) +
    ggtitle("b        ")
  
  layout4 <- rbind(c(1, 2),
                   c(3, 4))
  outname <- paste("../output/4paper/Fig4.png",sep="")
  png(filename=outname, width=3000, height=4500,res=500)
  grid.arrange(plot_fig4[[1]], plot_fig4[[4]], plot_fig4[[2]], plot_fig4[[3]],
               layout_matrix = layout4)
  dev.off()  

  plotdata0 <- plotdata1 %>% filter(Model %in% c("CGE1") & CLP %in% c("Mitigation","Mitigation_NoCCS","Mitigation_NoNuc")) %>% select(-SCENARIO)
  plot_fig4[[5]] <- 
    ggplot() + 
    geom_line(data=plotdata0,aes(x=Y, y = Value , color=CLP,group=interaction(CLP)),size=1.0) +
    geom_point(data=plotdata0,aes(x=Y, y = Value , color=CLP,group=interaction(CLP),shape=CLP),size=3.0,fill="white") +
    MyThemeLine+ scale_color_manual(values=linepalette) +
    scale_shape_manual(values=c(1,2,3)) +
    #    scale_shape_manual(values=ShapeVecSce) +
    xlab("Year") + ylab(paste0("Equivalent variation (",gdplossname2[,3],")") )  +
    annotate("segment",x=2005,xend=2050,y=0,yend=0,linetype="dashed",color="grey")+ 
    theme(legend.position="bottom", text=element_text(size=12), legend.title=element_blank(),
          axis.text.x=element_text(vjust=0.9, hjust=1, size = 12)) +
    guides(shape=guide_legend(ncol=1),color=guide_legend(ncol=1)) 
  outname <- paste("../output/4paper/Fig4_ref.png",sep="")
  ggsave(plot_fig4[[5]], file=outname, dpi = 250, width=5, height=5,limitsize=FALSE)
  
# GDP and popluation
  plot_SIfig10 <- as.list(2)
  silist1 <- c(34,33)
  alphabet <- c("a","b","c","d")
  plotdata1 <- allmodel %>% filter(Variable %in% c("GDP_MER","POP") & Model == "CGE1" & CLP=="Baseline") 
  for(i in 1:length(silist1)){
    plot_SIfig10[[i]] <- 
    ggplot() + 
    geom_line(data=plotdata1[plotdata1$Variable==varlist[silist1[i],1],c(-100)],aes(x=Y, y = Value),color="blue",size=1.0) +
    geom_point(data=plotdata1[plotdata1$Variable==varlist[silist1[i],1],c(-100)],aes(x=Y, y = Value),color="blue",size=3.0,fill="white") +
    MyThemeLine +
    xlab("Year") + ylab(paste0(varlist[silist1[i],2]," (",varlist[silist1[i],3],")") ) +
    annotate("segment",x=2005,xend=2050,y=0,yend=0,linetype="dashed",color="grey")+ 
    theme(legend.position="bottom", text=element_text(size=12), legend.title=element_blank(),
          axis.text.x=element_text(vjust=0.9, hjust=1, size = 12)) +
    ggtitle(paste0(alphabet[i],"      "))
  }
  
  plot_SIfig10[[2]] <- 
    ggplot() + 
    geom_line(data=plotdata1[plotdata1$Variable=="GDP_MER",c(-100)],aes(x=Y, y = Value),color="blue",size=1.0) +
    geom_point(data=plotdata1[plotdata1$Variable=="GDP_MER",c(-100)],aes(x=Y, y = Value),color="blue",size=3.0,fill="white") +
    MyThemeLine +
    xlab("Year") + ylab(paste0("GDP (million 2005US$/year)") ) +
    annotate("segment",x=2005,xend=2050,y=0,yend=0,linetype="dashed",color="grey")+ 
    theme(legend.position="bottom", text=element_text(size=12), legend.title=element_blank(),
          axis.text.x=element_text(vjust=0.9, hjust=1, size = 12)) +
    ggtitle("b      ")
  
  
  layout5 <- rbind(c(1, 2),
                   c(1, 2))
  #layout5 <- c(1, 2)
  outname <- paste("../output/4paper/SI_Fig4.png",sep="")
  png(filename=outname, width=4000, height=2000,res=500)
  grid.arrange(plot_SIfig10[[1]], plot_SIfig10[[2]],
               layout_matrix = layout5)
  dev.off()  
  
  #SI
#  indicator_names <- data.frame(c('CNS','Out_I_S','Out_Oth_Ind','Out_Ser'),c("Consumption","Iron and Steel","Other manufacturing","Service"))
  indicator_names <- data.frame(c('CNS','Out_I_S','Out_Oth_Ind','Out_Ser'),c("a                             ","b                             ","c                             ","d                            "))
  names(indicator_names) <- c("VEMF","Indicator")
  Mitbase_ratio <- rgdx.param('../output/temp.gdx','Mitbase_ratio') %>% filter(VEMF %in% c("CNS","Out_I_S","Out_Ser","Out_Oth_Ind"))%>% mutate(MODEL="CGE") %>% inner_join(indicator_names,by="VEMF")
  Mitbase_ratio$Y <- as.numeric(levels(Mitbase_ratio$Y))[Mitbase_ratio$Y]
  levels(Mitbase_ratio$VEMF) <- c("a", "b", "c","d")
  g <- ggplot()+
    geom_line(data=Mitbase_ratio,aes(x=Y, y = Mitbase_ratio*100,group=SMODEL,color=SMODEL),size=1.5) +
    geom_point(data=Mitbase_ratio,aes(x=Y, y = Mitbase_ratio*100,group=SMODEL,color=SMODEL),shape=21,size=2.0,fill="white") +
    scale_color_manual(values=pastelpal) +
    ylab(expression("Change ratios of each sector in mitigation \n compared to baseline scenario (%)")) + xlab("Year") +
#    scale_shape_manual(values=c(1,18,19)) +
    MyThemeLine 
  g2 <- g + facet_wrap(~Indicator,nrow=2)  +
  theme(legend.title=element_blank()) +
  annotate("segment",x=2005,xend=2050,y=0,yend=0,linetype="dashed",color="grey") + theme(legend.position='bottom') +
  theme(plot.margin = unit(c(0.5,0.5,0.5,0.5),"cm"))
  outname <- paste("../output/4paper/SIFig5.png",sep="")
  ggsave(g2, file=outname, dpi = 250, width=5, height=5,limitsize=FALSE)
  
## Iteration convergence
  iterationplot <- read.table("../data/iterationvar.txt", sep="\t",header=F, stringsAsFactors=F)
  iterationmap <- read.table("../data/iteration.map", sep="\t",header=T, stringsAsFactors=F)
  CLPlist <- as.vector(unique(allmodel.2$CLP))
  for(i in 1:length(CLPlist)){
    for(j in 1:2){
      yset <- c(2030,2050)
      allmodel.3 <- inner_join(allmodel.2,iterationmap,by="Model") %>% filter(Variable %in% as.vector(iterationplot$V1) & CLP==CLPlist[i]) 
      plot.multi.0 <- ggplot() + 
      geom_line(data=allmodel.3[allmodel.3$Y==yset[j],c(-2)],aes(x=Iteration, y = Value , color=MODEL,group=MODEL),stat="identity") +
      geom_point(data=allmodel.3[allmodel.3$Y==yset[j],c(-2)],aes(x=Iteration, y = Value , color=MODEL,group=MODEL,shape=MODEL),size=3.0) +
      MyThemeLine+ 
      scale_color_manual(values=linepalette)+ 
      #  scale_color_manual(values=linepalette,labels = c(CGE1="CGE Stand-alone",CGE2="Integrated model",Enduse1="Enduse Stand-alone")) + 
      #    scale_shape_manual(values=c(1,2,3,23)) +
      #  labels = c(CGE1="CGE Stand-alone",CGE2="Integrated model",Enduse1="Enduse Stand-alone")
      
      xlab("Iteration") + ylab("Unit depends on variables") +
    annotate("segment",x=0,xend=5,y=0,yend=0,linetype="dashed",color="grey")+ 
      theme(legend.title=element_blank()) 
    plot.1 <- plot.multi.0 + facet_wrap(~ Variable,scales="free")
    outname <- paste("../output/iteration/Iteration_",yset[j],"_",CLPlist[i],".png",sep="")
    ggsave(plot.1, file=outname, dpi = 250, width=16, height=12,limitsize=FALSE)
  }
  }
  
## Figure XXX
  varatio <- rgdx.param('../output/temp.gdx','GDP_share')
  Capprod <- rgdx.param('../output/temp.gdx','CapitaleffRel')
  capprodeff <- rgdx.param('../output/temp.gdx','VAChangeCapEff1')
  activitytop95 <- read.table("../data/activity_top95.txt", header=T, sep="\t",stringsAsFactors=F)
  modelmap <- read.table("../data/modelname.txt", header=T, sep="\t",stringsAsFactors=F)
  poweractivity <- read.table("../data/poweractivity.txt", header=T, sep="\t",stringsAsFactors=F)
    
  Capprod.1 <- Capprod %>% filter(SMODEL %in% c("CGE1","CGE2") & Y==2050) %>% inner_join(activitytop95,by="A")
  varatio.1 <- varatio %>% filter(SMODEL %in% c("CGE1","CGE2") & Y==2050) %>% inner_join(activitytop95,by="A")
  capprodeff.1 <- capprodeff %>% filter(SMODEL %in% c("CGE1","CGE2") & Y==2050) %>% inner_join(activitytop95,by="A") %>% inner_join(modelmap,by="SMODEL")
  capprodeff.tot <- capprodeff %>% filter(SMODEL %in% c("CGE1","CGE2") & Y==2050 & A=="TOT") %>% inner_join(modelmap,by="SMODEL")
  power.1 <- varatio %>% filter(SMODEL %in% c("CGE1","CGE2") & Y==2050) %>% inner_join(poweractivity,by="A") %>% inner_join(modelmap,by="SMODEL")
  vaplot <- as.list(1,2,3,4)
  
  vaplot[[1]] <- ggplot() +
    geom_line(data=Capprod.1,aes(x=Sector,y=CapitaleffRel,color=SMODEL,group=SMODEL),size=1.5) +
    geom_point(data=Capprod.1,aes(x=Sector,y=CapitaleffRel,color=SMODEL),shape=21,size=2.0,fill="white") +
    MyThemeLine+ 
    scale_color_manual(values=linepalette,labels = c(CGE1="CGE Stand-alone",CGE2="Integrated model")) +
    xlab("Sectors") + ylab("relative ratio of capital efficiency (-)") +
    ylim(0.85,1.2) +
    annotate("segment",x=0.5,xend=10.5,y=1,yend=1,linetype="dashed",color="grey")+ 
    theme(legend.title=element_blank()) 
  outname <- paste0("../output/valueadded/CapPro",".png")
  ggsave(vaplot[[1]], file=outname, dpi = 250, width=5, height=5,limitsize=FALSE)
  
  vaplot[[2]] <- ggplot() +
    geom_bar(data=varatio.1[varatio.1$SCENARIO=="Baseline" & varatio.1$SMODEL=="CGE1",c(-20)],aes(x=Sector,y=GDP_share,fill=SMODEL),stat="identity",fill=linepalette[1]) +
#    geom_point(data=Capprod.1,aes(x=A,y=CapitaleffRel,color=SMODEL),shape=21,size=2.0,fill="white") +
    MyThemeLine+ 
    xlab("Sectors") + ylab("Value added share in GDP (-)") +
    annotate("segment",x=0.5,xend=10.5,y=0,yend=0,linetype="dashed",color="grey")+ 
    theme(legend.title=element_blank()) 
  outname <- paste0("../output/valueadded/share",".png")
  ggsave(vaplot[[2]], file=outname, dpi = 250, width=4, height=4,limitsize=FALSE)
  
  ylabel <- expression("Contribution of capital productivity change to value added changes (% of GDP)")
  vaplot[[3]] <- ggplot() +
    geom_bar(data=capprodeff.1,aes(x=Model,y=VAChangeCapEff1*100,fill=Sector),stat="identity") +
    geom_point(data=capprodeff.tot,aes(x=Model,y=VAChangeCapEff1*100),shape=21,size=5.0,fill="white") +
    MyThemeLine+ 
    scale_fill_manual(values=linepalette)+ 
    xlab("Models") + ylab(ylabel) +
    annotate("segment",x=0.5,xend=2.5,y=0,yend=0,linetype="dashed",color="grey")+ 
    theme(legend.title=element_blank(),
          axis.text.y = element_text(size = 8)
    ) +
    coord_flip()
  outname <- paste0("../output/valueadded/CapProContribution",".png")
  ggsave(vaplot[[3]], file=outname, dpi = 250, width=7, height=4,limitsize=FALSE)
  
  vaplot[[4]] <- ggplot() +
    geom_bar(data=power.1,aes(x=interaction(SCENARIO,Model),y=GDP_share*100,fill=reorder(Power,-Order)),stat="identity") +
    MyThemeLine+         theme(legend.title=element_blank()) +
    scale_fill_manual(values=linepalette)+ 
    xlab("") + ylab(paste0(expression("Electricity sector's value-added \n share in GDP (%)"))) +
    annotate("segment",x=0.5,xend=4.5,y=0,yend=0,linetype="dashed",color="grey")
  outname <- paste0("../output/valueadded/va_power",".png")
  ggsave(vaplot[[4]], file=outname, dpi = 400, width=5.75, height=4.5,limitsize=FALSE)
  
  
#Elsticity of material and manetary translation variations
  plotdata.ev <- allmodel.2 %>% filter (Variable==gdplossname[,1])  %>% filter(Model %in% c("CGE1","CGE2E0p5","CGE2E2p0","CGE2","Enduse1","Enduse2","Enduse2E0p5","Enduse2E2p0") & CLP==c("Mitigation"))
#  plotdata.ev <- allmodel.2 %>% filter (Variable==gdplossname[,1])  %>% filter(Model %in% c("CGE1","CGE2E0p5","CGE2E2p0","CGE2") & CLP==c("Mitigation"))
  plotdata.ev <- na.omit(plotdata.ev)%>% select(-SCENARIO)
  plot_fig_ev <- 
    ggplot() + 
    geom_line(data=plotdata.ev,aes(x=Y, y = Value , color=Model,group=interaction(Model)),size=1.0) +
    geom_point(data=plotdata.ev,aes(x=Y, y = Value , color=Model,group=interaction(Model)),shape=21,size=3.0,fill="white") +
    MyThemeLine+
#    scale_shape_manual(values=c(1,2,3,4,23)) +
    scale_color_manual(values=linepalette,labels = c(CGE1="CGE Stand-alone",CGE2="Integrated model",CGE2E0p5="Integrated model (ela=0.5)",CGE2E2p0="Integrated model (ela=2.0)",
                                                     Enduse1="Enduse Stand-alone",Enduse2="Enduse second",Enduse2E0p5="Enduse second (ela=0.5)",Enduse2E2p0="Enduse second (ela=2.0)")) + 
    #    scale_shape_manual(values=ShapeVecSce) +
    xlab("Year") + ylab(paste0("Policy cost (",gdplossname[,3],")") ) +
    annotate("segment",x=2005,xend=2050,y=0,yend=0,linetype="dashed",color="grey")+ 
    theme(legend.position="bottom", text=element_text(size=12), legend.title=element_blank(),
          axis.text.x=element_text(vjust=0.9, hjust=1, size = 12)) +
    guides(shape=guide_legend(ncol=2),color=guide_legend(ncol=2)) 
  outname <- paste0("../output/4paper/elasticityvariant",".png")
  ggsave(plot_fig_ev, file=outname, dpi = 400, width=4, height=5,limitsize=FALSE)
  
#---- For paper figure end 
  
